Controlling the Growth of Constraints in Hyperbolic Evolution Systems 
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Motivated by the need to control the exponential growth of constraint violations in numerical 
solutions of the Einstein evolution equations, two methods are studied here for controlling this 
growth in general hyperbolic evolution systems. The first method adjusts the evolution equations 
dynamically, by adding multiples of the constraints, in a way designed to minimize this growth. 
The second method imposes special constraint preserving boundary conditions on the incoming 
components of the dynamical fields. The efficacy of these methods is tested by using them to control 
the growth of constraints in fully dynamical 3D numerical solutions of a particular representation of 
the Maxwell equations that is subject to constraint violations. The constraint preserving boundary 
conditions are found to be much more effective than active constraint control in the case of this 
Maxwell system. 



I. INTRODUCTION 

Perhaps the most critical problem facing the numerical 
relativity community today is the exponential growth of 
constraints in evolutions of black hole spacetimes. The 
evolution equations guarantee that constraints that are 
satisfied exactly on a spacelike surface will be satisfied 
throughout the domain of dependence of that surface. 
However, this result does not guarantee that small initial 
violations of the constraints will remain small, or that 
constraint violations will not be injected into the com- 
putational domain through timelike boundaries. Experi- 
ence has shown that constraint violations tend to grow 
exponentially in the numerical evolution of black hole 
spacetimes (e.g., 0, H 0)- These constraint violating 
instabilities have been shown to be numerically conver- 
gent and thus represent unstable solutions to the partial 
differential equations. At present these instabilities are 
the limiting factor that prevents these numerical simula- 
tions from running for the needed time with the required 
accuracy. 

A variety of approaches have been explored in a num- 
ber of attempts to control the growth of the constraints: 

1) Fully constrained evolution, in which the constraint 
equations are re-solved periodically (e.g. at each time 
step) have been used with great success in spherically 
symmetric and axisymmetric problems 0, H, El B HI 
llCj . These methods have not gained widespread use in 
3D simulations, however, due in part to the high cost of 
solving the elliptic constraint equations. Difficult ques- 
tions also remain unresolved for this method about the 
proper boundary conditions to impose on the constraint 
equations at black hole excision boundaries. With the 
development of more efficient elliptic solvers and the ab- 
sence of a better alternative however, fully constrained 
methods are starting to be developed and tested in 3D 
now as well E3. 

2) Auxiliary dynamical fields have been introduced into 
the system whose evolution equations are designed to 
drive the constraints toward zero [blj . This technique 



has the disadvantage that it requires the size of the dy- 
namical system to be significantly expanded. It has not 
been tested extensively, but the first numerical results 
were not uniformly successful |l5l llq . 

3) More sophisticated boundary conditions have been 
introduced whose purpose is to control the influx of con- 
straint violation throu gh t he timelike boundaries of the 
computational domain [H mm^SIHE!!!!!! 
[2(J. This approach seems very promising, although the 
current methods may not be fully compatible with the 
physical requirement that waves pass through the bound- 
aries without reflection. Further these boundary condi- 
tion methods may not completely solve the constraint 
violating instability problem in systems like the Ein- 
stein evolution equations, where constraint violations are 
driven both by bulk and by boundary terms in the equa- 
tions. But this technique can (as we will demonstrate 
below) significantly improve the influx of constraint vio- 
lations through the timelike boundaries of the computa- 
tional domain. 

4) Dynamically changing the evolution equations, 
through the addition of terms proportional to the con- 
straints, has been proposed as a way to minimize con- 
straint growth. This method (developed by Manuel 
Tiglio and his collaborators [53, ) has had some suc- 
cess in controlling the growth of constraints in simple nu- 
merical solutions of the Einstein evolution equations. We 
find that this technique when used in combination with 
standard boundary conditions is not effective however in 
controlling the influx of constraint violations through the 
boundaries of the computational domain in fully dynam- 
ical situations. 

In this paper we explore in some detail two of these 
methods for controlling the growth of constraint viola- 
tions in hyperbolic evolution systems. First, we develop 
a refined version of the dynamical constraint control 
method being used by Tiglio and collaborators [27L |2S| . 
In particular we introduce a more natural norm on the 
constraints, which has the property that its evolution 
can be predicted numerically with greater accuracy. We 
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expect that dynamical constraint control based on this 
new constraint norm should be more stable and robust 
than the current method. And second, we explore the 
use of constraint preserving boundary conditions. In this 
method (exp lored previously by Calabrese and collabora- 
tors [TflL I2{t|) the constraints are decomposed into char- 
acteristic ingoing and outgoing fields of the constraint 
evolution equations. Setting the incoming components of 
the constraint fields to zero provides boundary conditions 
for the incoming parts of the dynamical fields of the prin- 
cipal evolution system. We test both of these methods 
by applying them to a non-trivial hyperbolic evolution 
system (a p articular representation of the Maxwell sys- 
tem [23, that is analogous to, but much simpler than, 
the Einstein evolution system. Our tests — using fully dy- 
namical time dependent solutions on domains with open 
boundaries — reveal that the constraint preserving bound- 
ary conditions are much more effective than active con- 
straint control for this Maxwell system. Some features 
of this system are rather special, and it is possible that 
in more generic systems (like the Einstein equations) the 
active constraint control method may be complimentary 
to the constraint preserving boundary condition method. 



We define and review in Sec. UTI t he particular form 
of the Maxwell evolution system |29|, |30| that we use to 
illustrate and test the constraint control methods stud- 
ied here. We refer to this system as the "fat" Maxwell 
system since it replaces the usual representation of the 
Maxwell system, which has six independent field com- 
ponents, with a representation having twelve. We also 
present in Sec. [D] the decomposition of the dynamical 
fields used in this fat Maxwell system into characteris- 
tic parts. In Sec. lIIII we develop the equations needed to 
perform active constraint control, in particular on the fat 
Maxwell evolution system. We determine the constraint 
evolution equations for this system, and derive an im- 
proved norm on the constraint fields. We show how the 
evolution of this new constraint norm should generically 
be more accurately determined (and hence should pro- 
vide better constraint control) in numerical solutions. In 
Sec. lIVI we develop the particular form of constraint pre- 
serving boundary conditions studied here. We present 
the decomposition of the constraint fields into character- 
istic parts, and show how these can be used to provide 
boundary conditions for the principal evolution system. 
Finally in Sec. we use these methods to control the 
growth of constraints in fully dynamical 3D numerical 
evolutions of the fat Maxwell evolution system. We note 
that both the active constraint control mechanism and 
the constraint preserving boundary conditions developed 
here are applicable to rather general hyperbolic evolution 
systems. We focus our discussion on the fat Maxwell sys- 
tem in order to make the analysis less abstract, and to 
provide a simple system on which to perform numerical 
experiments. 



II. FAT MAXWELL EVOLUTION SYSTEM 

Our primary interest here is to understand how to con- 
trol the growth of constraints in hyperbolic evolution sys- 
tems. We will focus our attention on quasi-linear systems 
of the form, 

d t u a +A ka l3 d k u> 3 =F Q , (1) 

where u a represents the dynamical fields, and A ka p and 
F a may depend on u a but not its derivatives. We assume 
that the evolution system described in Eq. Q is also 
subject to a set of constraints, c A = 0, which we assume 
have the general form 

c A = K Ak a d k u a + L A , (2) 

where K Ak p and L A may depend on the dynamical fields 
u a but not their derivatives. We assume that these con- 
straints are preserved as a consequence of the evolution 
equations. In particular we assume that the constraints 
satisfy an evolution equation of the form 

d t c B + A kB D d k c D = F B D c D , (3) 

where A kB o may depend on the dynamical fields u a , 
while F B d may depend on u a and its spatial deriva- 
tives d k u a . When this constraint evolution system is hy- 
perbolic the constraints will remain satisfied within the 
domain of dependence of the initial surface if they are 
satisfied initially. We note that multiples of the con- 
straints of the form given in Eq. J2J may be added to 
the principal evolution system Eq. without changing 
the physical (constraint satisfying) solutions of the sys- 
tem or the basic structure of Eq. . Systems with this 
general form include most of the evolution equations of 
interest in mathematical physics, including for example 
the Einstein evolution equations, the Maxwell equations, 
the incompressible fluid equations, etc. 

In order to explore and test some of the ideas for con- 
trolling the growth of constraints in these hyperbolic evo- 
lution systems, we adopt a simple example system on 
which to perform our analysis and to carry out numeri- 
cal tests. We have chosen to use a form of the vacuum 
Maxwell evolution equations (introduced independently 
by Kidder [2^| and Reula [3(|) that fits nicely into this 
framework, and that admits constraint violations if noth- 
ing is done to control them. The dynamical fields in this 
formulation are a co-vector that represents the electric 
field Ei, and a second rank tensor Z?y that represents 
the gradient of the spatial parts of the vector potential 
(i.e. Dij = diAj, although we impose the relationship 
between and the vector potential only implicitly as 
a constraint on this system). We refer to this as the 
fat Maxwell system, since the usual representation of 
the Maxwell equations with six dynamical field compo- 
nents is replaced with this larger representation that has 
twelve. The evolution equations for this system are, 

d t Ei = g ab d a (D lb ~D u ), (4) 
d t Dij = -diEj - didj(j>, (5) 
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where g ab is the Euclidean metric with inverse g ab , and d a 
is the covariant derivative compatible with this metric, 
(i.e. just partial derivatives in Cartesian coordinates). 
The scalar potential </> is a gauge quantity assumed here 
to be a given function of space and time. This system has 
the same general form as Eq. Q with u a = {Ei, Dij}. In 
order to represent the vacuum (i.e. charge and current 
free) Maxwell system these equations are also subject to 
the constraints, C = Ctj k = 0, where 



C = g ab d a E b , 



(6) 
(7) 



These constraints have the same general form as those 
described in Eq. (J2J with c A — {C,Cijk}- The second of 
these constraints is the integrability condition that guar- 
antees that is the gradient of a vector potential. As 
mentioned above we are free to add multiples of the con- 
straints to the evolution system: 

dtEi = g ab d a (D lb -D M ) + lig ab C iab , (8) 
d t Dij = -diEj - didj<j> + ~fygijC, (9) 

where 71 and 72 are constants. The addition of these 
constraint terms leaves the physical (constraint preserv- 
ing) solutions to the system unchanged, and also leaves 
the system with the same basic structure as Eq. Q . 

For hyperbolic evolution systems, such as those in 
Eq. (JTJ , it is often quite useful to decompose the dynam- 
ical fields u a into characteristic fields. These character- 
istic fields are defined with respect to a spatial direction 
at each point, represented here by the unit normal co- 
vector field rife. Given a direction field n k we define the 
eigenvectors e a a of the characteristic matrix A ka p: 



(10) 



where denotes the eigenvalue (also called the char- 
acteristic speed). The index a labels the various eigen- 
vectors and eigenvalues, and there is no summation over 
a in Eq. (|10|l . Since we are interested in hyperbolic evo- 
lution systems, the space of eigenvectors will have the 
same dimension as the space of dynamical fields, and the 
matrix e a p will be invertible. Given these characteristic 
eigenvectors it is often useful to re-express the dynam- 
ical fields in terms of this eigenvector basis. Thus we 
define the characteristic fields u a (or the characteristic 
projection of the dynamical fields) as 



(11) 



It is straightforward to show that the evolution of the 
characteristic fields is determined by 

d t u & +v (&) n k d k u & = -e & a P n k A ka pd n u? + e & a F a 

+ (d t e & a +v {&) n k d k e & a )u a ,(12) 

where the projection operator orthogonal to is defined 
by Pij — gij — riirij, and spatial indices are raised and 
lowered with g tJ and gij respectively. 



The characteristic fields for the fat Maxwell evolu- 
tion system are a collection of fields of the form u a = 



{Z 1 ,Zf,Zf i ,Ul ± ,U 2± },where 
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(17) 



The characteristic fields Z 1 , Zj and Zf 4 have character- 
istic speed 0; the fields have speeds ±1, and the 
fields U have speeds ±-^/7i72- All the characteristic 
speeds are real, and the characteristic fields are linearly 
independent (and depend continuously on the unit vector 
n k ) whenever 7172 > 0. Consequently the fat Maxwell 
system is strongly hyperbolic when 7172 > 0. We also 
find that the fat Maxwell evolution system is symmet- 
ric hyperbolic when the parameters 71 and 72 satisfy the 
more restrictive conditions, < 71 < 4, and | < 72. 

We note that the characteristic eigenvectors e a a for the 
fat Maxwell system depend only on the spatial metric g^ 
and the normal vector n,. It follows that the last term 
on the right side of Eq. I|12|) does not depend on any 
derivatives of u a at all. Thus the right side of Eq. i|12|) 
depends only on the transverse (to rii) derivatives of u a : 



d t u a + v (&) n k d k u a = G a {u l3 ,P k n d k u p ). 



(18) 



This important feature of the characteristic evolution 
equations is satisfied by many systems of interest to us, 
including the Einstein evolution system. 

It is also useful to have the inverse transformation 
u a = e a aU a , where e a a is the inverse of e a a . For the 
fat Maxwell system this inverse transformation has the 
form: 



E 
D 



i(C// + + Ut) + in 4 ([/ 2+ -U 2 -), 

'Z 1 72-1 



(19) 
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72 



V7i72 

:(U 2+ - 



V 1 



2+ 



Z 2 rij 



VT02 

+ \ ni [U] + - Uj~ - (71 - 2)Z 2 } + Zfj. (20) 



The characteristic decomposition of the dynamical fields 
is essential for fixing boundary conditions. We will return 
to a more complete discussion of boundary conditions in 
Sec. El 



III. ACTIVE CONSTRAINT CONTROL 

Unless the constraint evolution system Eq. (PJ) is hy- 
perbolic, it will not guarantee that the constraints re- 
main satisfied (within the domain of dependence of an 
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initial surface) even if they are satisfied initially. Thus 
the constraint evolution system must be hyperbolic in 
any self-consistent and physically reasonable system of 
constrained evolution equations. We assume that any 
system considered here has a hyperbolic constraint evo- 
lution system. We also assume that the constraint evolu- 
tion system satisfies the somewhat stronger condition of 
symmetric hyperbolicity: In particular we assume that 
there exists a symmetric, positive-definite tensor Sab on 
the space of constraint fields which symmetrizes the char- 
acteristic matrices of the constraint system, 

S AC A kC B = A k AB = A k BA , (21) 

for all k. When such a symmetrizer exists, we can define 
a natural norm on the constraints: The constraint energy 
£ and its associated current £ k are defined by 

£ = S AB c A c B , (22) 
£ k = A k AB c A c B . (23) 

This constraint energy can be used to define a norm (£ ) 
on the constraints, 

(£) = J £d 3 x, (24) 

since (£) = if and only if all the constraints are satisfied 
at each point. It is straightforward to determine the time 
evolution of £ using the constraint evolution equations for 
any symmetric hyperbolic constraint evolution system: 

d t £ + d k £ k = £ AB c A c B . (25) 

The quantities £ k and £ ab may depend on the dynamical 
fields u a and their spatial derivatives d k u a (but not on 
higher spatial derivatives of u a ). 

In an evolution system Eq. Q that includes parame- 
ters 7 a multiplying constraint terms, such as the system 
defined by Eqs. JSJl and (JSJ, the associated constraint evo- 
lution system Eq. J3J and the constraint energy system 
Eq. H25fl will also include terms that depend linearly on 
these parameters. Integrating Eq. I|25[) over the spatial 
slice for such a system, we get an expression for the time 
evolution of the constraint norm which has the general 
form, 

d t (£) = Q + iaR a , (26) 

where Q and R a are integrals of quantities that depend 
on the dynamical fields and their first spatial derivatives. 
The basic idea of active constraint control then is to ad- 
just the parameters j a that appear in Eq. I|26|l to control 
the evolution of the constraint norm (£). For example 
the growth of (£ ) might be prevented by making the right 
side of Eq. (|2(j[) non-positive at the beginning of each time 
step in the numerical evolution. This control mechanism 
is a special case of the constraint control method devel- 
oped by Tiglio and his collaborators I27 L |2 8l . I t differs 
from Tiglio's particular implementation |27l |28| in that 



the quantities Q and R a in our expression do not de- 
pend on second derivatives of the dynamical fields. Since 
these higher derivatives are more difficult to evaluate ac- 
curately in a numerical simulation, we expect that our 
constraint control mechanism will be more stable and ro- 
bust. 

The constraints associated with the vacuum fat 
Maxwell system introduced in Sec. |n] satisfy the follow- 
ing evolution system as a consequence of Eqs. (fSJ and 
©, 

d t C = 7 i9 ii 5 a6 ^C i „ 6 , (27) 
d t Cijk = ^(fjjkdiC - gikdjC). (28) 

This system has the same general form as Eq. @ with 
c A = {C,Cijk}- In order to define a constraint energy, 
we need this constraint evolution system to be symmet- 
ric hyperbolic. The most general symmetrizer for this 
system (that can be constructed from the spatial metric 
g ab ) is given by 

dS 2 = S A Bdc A dc B 

= xig ij s/ ah dCi j de ab + X 2g ia g jb dC l] dC ab 
+ X 3g ia g lb dC M dC [ah] + 2 X3 ^dC 2 , (29) 

where 

dC t] = g lc e cab dC ab3 , (30) 
dCa = h{5i5) + 5 a A-Uii9 ab )dCab, (31) 

and e 1 ^ is the totally antisymmetric tensor volume ele- 
ment. The parameters \ a must be positive \ a > 0, and 
7i72 must also be positive 7172 > to ensure that Sab 
is positive definite. We note that these conditions put no 
additional limits on the allowed ranges of the parameters: 
every strongly hyperbolic representation of the principal 
evolution system has a symmetric hyperbolic constraint 
evolution system. 

We now evaluate the various quantities that determine 
the evolution of the constraint energy, Eq. 1)25(1 . for the 
fat Maxwell system. We find: 

£ k = -4 l2X 3C^g ka C alJ , (32) 
£ab = 0. (33) 

Thus the expression for the time derivative of the con- 
straint energy becomes, 

dt£ = 4j2X3d k (CgVg ka C alJ ). (34) 

The right side of Eq. I|34|) is a divergence, so the inte- 
gral of this equation over a spatial surface results in an 
expression that involves only boundary integrals: 

d t {£) = 4-Y2X3</>Cg ij n k C kij d 2 x, (35) 

where n k is the outward directed unit normal to the 
boundary. Active constraint control for this system con- 
sists then of adjusting the sign of the parameter 72 to 
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force the constraint norm (£) to decrease with time 
whenever it gets unacceptable large. 

We note that the fat Maxwell system is rather degener- 
ate, since the right side of Eq. I|35l) contains only a surface 
integral. Thus constraint violation in the fat Maxwell 
system arises only from the influx of constraint violations 
through the timelike boundaries of the computational do- 
main. This property makes this system rather simpler 
than the Einstein evolution equations where constraint 
violation can be generated from bulk terms in the equa- 
tions as well. The simplicity of the fat Maxwell system 
allows us to study how best to control the influx of con- 
straint violations across boundaries in some detail, but it 
does not allow us to evaluate how effective these methods 
are for controlling violations that arise from bulk terms 
in the equations. 

IV. CONSTRAINT PRESERVING BOUNDARY 
CONDITIONS 

A standard boundary condition used for hyperbolic 
systems is the maximally dissipative condition, which 
we define here to be the condition that sets the incom- 
ing components of the dynamical fields to zero. (More 
generally the term maximally dissipative has been used 
to describe a larger class of boundary conditions that 
guarantee that a certain energy flux at the boundaries is 
strictly outgoing, e.g. see 13.). To impose such a condi- 
tion, we first decompose the dynamical fields into char- 
acteristic parts, as was done in Eq. I|llfl . and then set to 
zero at the boundary all those characteristic fields whose 
characteristic speeds are negative. Let EPa denote the 
projection operator that annihilates all the non-incoming 
characteristic fields: that is, let 

n & y = !f i 0T ^)<°< (36) 

p [ for V(&) > 0. v ' 

For a maximally dissipative boundary condition, we set 
n a mP — at the boundaries. We often use a varia- 
tion on this boundary condition, in which we set to zero 
the time derivatives of the incoming components of the 
characteristic fields: 

\T^d t J = 0. (37) 

For the case of the fat Maxwell system discussed in 
Sec. m these "freezing" boundary conditions reduce to 

d t Ul~ = d t U 2 ~ = 0, (38) 

where the incoming characteristic fields Xj\~ and U 2 ~ are 
defined in Eqs. I(16|) and (| 1 T|l . As we shall see in Sec. Ivl 
this "freezing" boundary condition does a poor job of 
preventing the influx of constraint violations through the 
boundaries. 

Calabrese, et al. 19] have proposed an alternate 
method for constructing boundary conditions that pre- 
vent the influx of constraint violations. Their method 



involves decomposing the constraint fields c A into char- 
acteristic parts: 

c A = e A B c B , (39) 

where e A a represents the eigenvectors of the character- 
istic matrix of the constraint evolution system, 

e A B n k A kB c = v (A) e A c, (40) 

and v,£, represents the corresponding eigenvalue (or 
characteristic speed). The idea is to impose what 
amounts to maximally dissipative boundary conditions 
on the constraint evolution equations: that is, we set 

Tl A § c B = 0, (41) 

where H A g is the projection operator that annihilates the 
non-incoming characteristic constraint fields. This con- 
dition must now be converted into a boundary condition 
on the dynamical fields of the principal evolution system 
u a . This is done through the equation that defines the 
constraints in terms of u a and its derivatives, Eq. J3J). In 
particular we solve Eq. I|41|) for the normal derivatives of 
the incoming characteristic fields, in terms of the outgo- 
ing characteristic fields and tangential derivatives of the 
incoming fields. When this is possible, Eq. (|41|l becomes 
a Neumann-like boundary condition on (some of) the in- 
coming characteristic fields. This boundary condition has 
the following general form 

n k d k J = H & [J, {5^ - n^)d k u\n^P k n d k u\ (42) 

We illustrate this procedure below more explicitly (and 
perhaps more clearly) for the specific case of the fat 
Maxwell system. 

The characteristic fields for the fat Maxwell constraint 
system are the collection of fields of the form c A = 



{Zf,Zf p U 3± }, where 

Zf = C [ik] n k , (43) 

4 = %)' ( 44 ) 

U 3± = -^^-C±n k g^C kl3 . (45) 



The fields Zf and Zfj have characteristic speed 0, while 
the fields U 3± have speeds ±^/7i72- The only incoming 
characteristic field is U 3 ~. So the constraint preserving 
boundary condition sets U 3 ~ = on the boundaries of 
the computational domain. Using the definition of U 3 ~ 
above, we see that this boundary condition is equivalent 
to the condition 

«V^--^C (46) 
7i 

on the boundaries. For a solution that satisfies the con- 
straint preserving boundary condition, Eq. 1(46(1 . the evo- 
lution of the constraint energy norm Eq. 1)35(1 becomes 

d t (£) = -A X3 Vl^^fc 2 d 2 x<0. (47) 
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Thus the constraint preserving boundary condition en- 
sures that the constraint norm does not grow. Quite 
generally the constraint preserving boundary conditions 
of this type will ensure that surface flux terms do not 
contribute to the growth of the constraint energy. 

In order to convert the constraint preserving bound- 
ary condition into an explicit condition on the dynamical 
fields, we must express the incoming constraint field U 3 ~ 
in terms of the characteristic fields u a . Using Eqs. I|13|l - 
(|T7|) and (gSJ) w e obtain 

u 3 - = ^ni[n k d k u 2 - - ^'a i (f/ 7 1+ + u)-)\ 

7i 

+ kP ij di [Uj+ - Uj- - (71 - 2)Zf] . (48) 

Setting U 3 ~ — - we obtain an expression for the normal 
derivative of U 2 ~: 

n k d k U 2 - = iP ij di(Uj + + Uj') 

i-^=P ij di [Uj + -Uj--(yx- 2)Z 2 } . (49) 
V7i72 

This has the form of a Neumann-like boundary condition 
on U 2 ~ , and has the same form as the general expression 
Eq. ftSj . 

The version of our code used to perform the numerical 
tests described in Sec. Ivl imposes boundary conditions on 
the time derivatives of the incoming characteristic fields. 
We therefore convert the Neumann-like boundary con- 
dition on U 2 ~ in Eq. I|49l) into a condition on its time 
derivative using the characteristic field evolution equa- 
tion, Eq. (|12[1 . We simply replace the normal derivative 
n k dkU 2 ~ that appears in Eq. Ijl2(l with the expression 
from Eq. (|49|l . Simplifying the results gives the following 
equation for the time derivative of U 2 ~ at the boundary, 

d t U 2 - = i^^-P^idiEj + didjcf) 
7i 

+2P^n k d l D m . (50) 

This condition together with the freezing boundary con- 
ditions dtU\~ on the remaining incoming characteris- 
tic fields constitute our version of constraint preserving 
boundary conditions on the fat Maxwell system. 

V. NUMERICAL RESULTS 

In this section we present numerical experiments that 
illustrate the effectiveness of the various constraint con- 
trol strategies discussed in this paper. All of these tests 
use the fat Maxwell evolution system, with a variety of 
topologies for the computational domain and with a va- 
riety of initial data for the dynamical fields. In Sec. IV Al 
we illustrate what happens when the equations are solved 
without any constraint control. These tests show that 
significant constraint violations (and in some cases con- 
straint violating instabilities) occur in dynamical solu- 
tions of the fat Maxwell system on computational do- 
mains with open boundaries. In Sec. IV Bl we study the 



use of the active constraint control mechanism described 
in Sec. 11111 Our tests show that this method is not nu- 
merically convergent, and is not very effective in con- 
trolling the growth of constraints in this system. And 
finally in Sec. IV CI we describe the results of using the 
constraint preserving boundary conditions described in 
Sec. IIVI This method is shown to be numerically con- 
vergent and quite effective in controlling the growth of 
constraints in the symmetric hyperbolic subset of the fat 
Maxwell system. 

All numerical computations presented here are per- 
formed using a pseudospectral collocation method. Our 
numerical methods are essentially the same as those we 
have applied to the evolution problem in full general rel- 
ativity 0, EE El OOI an d in studies of scalar fields in Kerr 
spacetime 32]. Given a system of partial differential 
equations 

d t u a (*,t) = .r>(x,i),d lU (x,i)], (51) 

where u a is a vector of dynamical fields, the solution 
u a (x, t) is expressed as a time-dependent linear combi- 
nation of N spatial basis functions 0fc(x): 

N-l 

i#(x,t) = J>g(t)&(x). (52) 

k=0 

Spatial derivatives are evaluated analytically using the 
known derivatives of the basis functions: 

AT-l 

c^(x,t) = ^<(i)c^(x). (53) 

k=0 

The coefficients u k (t) are chosen so that Eq. I|51|l is sat- 
isfied exactly at N c collocation points x< selected from 
the spatial domain. The values of the coefficients are 
obtained by the inverse transform 

w c -i 

= J2 u%(xi,t)</> k (xi)wi, (54) 

i=0 

where Wi are weights specific to the choice of basis func- 
tions and collocation points. It is straightforward to 
transform between the spectral coefficients u k (t) and the 
function values at the collocation points u^(xj,f) using 
Eqs. (|52|l and i|54[l . The partial differential equations, 
Eq. I|51|) . are now rewritten using Eqs. (|52l) -il54 [) as a set 
of ordinary differential equations for the function values 
at the collocation points, 

d t u%(x i ,t)=g?[u N (x d ,t)], (55) 

where Qf depends on u%{n.j, t) for all j. This system of 
ordinary differential equations, Eq. (|55l) . is integrated in 
time using a fourth-order Runge-Kutta method. Bound- 
ary conditions are incorporated into the right-hand side 
of Eq. ijSlH) using the technique of Bj0rhus (33|. The 
time step is typically chosen to be half the distance be- 
tween the closest collocation points, which ensures that 
the Courant condition is satisfied. 



7 



In order to provide a quantitative measure of conver- 
gence and the amount of constraint violation of our nu- 
merical solutions, we have defined several norms on the 
constraints c A and the dynamical fields u a . We have al- 
ready defined the constraint energy (£ ) in Eq. I|24l) , which 
provides a norm on the constraint fields. In computing 
(£) for these numerical studies we fix \i = Xi — X3 — !• 
We also define norms on the dynamical fields themselves: 



(E t E l 



D tj D I: >)d 3 x, 



\u\\ 2 LOO = max(i?i£ , '' + D i3 D ij ) 



(56) 
(57) 



We compute the volume integrals in these norms, e.g. 
in Eq. iffi^jl or (JSSJ, exactly using the appropriate form 
of Gaussian quadrature, and the maximum in Eq. i)57|) 
is taken over the appropriate set of collocation points 
at a particular instant of time. These norms are most 
useful for comparing solutions evaluated with different 
numerical resolutions. Thus we define 



INI 



J {6E.5E 1 + ) d 3 x, (58) 

(59) 



where Su a = {5Ei,8Dij} is the difference between u a 
at a given resolution and u a at the best (highest) res- 
olution we computed. Differences between quantities 
at different resolutions are computed by evaluating and 
then subtracting the spectral series for each resolution 
at the points on the finer grid. In order to provide 
meaningful scales for these normed quantities we typ- 
ically plot dimensionless ratios of expressions such as 
INIia/llubertHia and INHoo/IKestlHoo. In the case 
of the constraint energy we typically plot (£)/||9u|| 2 , 
where 



\\du\f 



(60) 



is a norm on the gradients of the fields. We are interested 
in seeing how these ratios behave as the resolution of 
the numerical solution is increased: Order unity values 



of these ratios, \\Su\ 



|u b est|| 2 or {£}/\\du\\ 2 , indicate 



a complete lack of numerical convergence or solutions 
that are dominated by constraint violations, respectively. 
Values of these ratios of order 10~ 34 correspond to double 
precision roundoff error. 



A. No Constraint Control 

In this section we illustrate the results of finding nu- 
merical solutions to the fat Maxwell evolution system 
Eqs. using no constraint control at all. We exam- 

ine three separate cases: First we look at evolutions on a 
computational domain with topology T 3 , a 3-torus. The 
differential equations governing the fat Maxwell system 
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FIG. 1: Convergence test for fat Maxwell on T . Shown are 
norms of differences between solutions at different resolutions: 
solid curves use the L 2 norms while the dotted curves use the 
L°° norms. 




20 30 
Time 



FIG. 2: Constraint violation for fat Maxwell on T . Shown is 
the constraint energy (£ ) divided by the norm of the deriva- 
tives of the fundamental variables. 



allow no constraint growth on domains without bound- 
aries. So this first test is to verify that our code accu- 
rately reproduces "exact" constraint conservation in this 
case. Next we examine the evolution of a representation 
of the static Coulomb solution on a computational do- 
main with topology S 2 x R, a spherical shell. Finally we 
study a highly dynamical solution on a computational 
domain with topology S 2 x R using freezing boundary 
conditions and no constraint control. 

The evolution of the constraint energy norm (£) for 
the fat Maxwell system is driven entirely by a boundary 
term, Eq. H35fl . Thus we expect the constraints to be 
satisfied exactly for evolutions on a computational do- 
main without boundary. To confirm that our numerical 
code correctly models this, we solve Eqs. |[5)-|[5J on a 
computational domain with topology T 3 , i.e. within a 
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3-torus. In particular we choose coordinates x, y, and 
z in the interval [0, 2n], and impose periodic boundary 
conditions. The basis functions used in our pseudospec- 
tral method are sines and cosines. We use initial data for 
this test in which the electric field is set to zero, Ei = 0, 
and each component of the vector potential is set to be 
a cylindrical Gaussian pulse: 

A x = Ay = A z = e - [{y - c y )2+{z - c * )2]/w \ (61) 

where the width of the pulse is set to w = 0.5, and the 
center is placed in the middle of the computational do- 
main, c y = c z — 3.14. The shape of this pulse is selected 
so that the value of the Gaussian falls below double pre- 
cision roundoff, 1CP 17 , at the periodicity "boundaries" of 
the domain, y = and y = 2ir etc. This ensures that 
these data are smooth on T 3 to roundoff accuracy. The 
initial data for Dij are set to the numerically determined 
values of diAj . We use the gauge choice <p = throughout 
this evolution. Because these initial data are effectively 
two-dimensional, we can place as few as two collocation 
points in the x direction for computational efficiency. 

Figurenshows a convergence plot for this case that was 
run with evolution parameter values 71 = I/72 = —0.1, 
and resolutions N y = N z = 10,20,30,40, and 50 col- 
location points. We see from Fig. 2] that the differences 
converge to zero as the resolution is increased. FigureElil- 
lustrates the amount of constraint violation in these runs. 
These curves, which increase approximately linearly with 
time, have magnitudes that are roughly proportional to 
the number of numerical operations performed multiplied 
by double precision roundoff error. Thus, the finer res- 
olutions have larger errors than the coarser ones, since 
the finer resolutions require a larger number of timesteps 
and a larger number of numerical operations at each step. 
As expected from Eq. l|35|l . we see that the constraints 
are satisfied essentially exactly when the domain has no 
boundaries. We have also computed evolutions for these 
initial data on T 3 using other values of the evolution 
parameters. In particular we have computed evolutions 
with 71 = I/72 =0.1, and also evolutions that switch 
back and forth between these positive and negative val- 
ues at each time step. In all of these cases, we find the 
evolutions to be convergent with roundoff level constraint 
violation. 

Next we turn our attention to solving the evolution 
equations on a computational domain with topology 
S 2 x R, i.e. within a spherical shell. For our basis func- 
tions we choose Chebyshev polynomials for the radial 
coordinate and spherical harmonics for the angular co- 
ordinates. Although our basis functions are written in 
(r, 9, ip) coordinates, our fundamental variables are the 
Cartesian components of Ei and Dij . To eliminate high- 
frequency numerical instabilities that sometimes develop 
during our simulations in S 2 x i?, we apply a filter to the 
right-hand sides of Eqs. i|55|) before and after incorporat- 
ing boundary conditions via the Bj0rhus algorithm. The 
filter consists of simply setting high-frequency spherical 
harmonic coefficients to zero: If £ ma x is the largest index 
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FIG. 3: Convergence test for fat Maxwell on S 2 x R with 
static point charge initial data. Shown are norms of differ- 
ences between solutions at different resolutions: solid curves 
use L 2 norms and dotted curves use L°° norms. 
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FIG. 4: Constraint violation for fat Maxwell on S 2 x R with 
static point charge initial data. Shown is the constraint en- 
ergy {£) divided by the norm of the derivatives of the funda- 
mental variables. 



used in the basis functions Yg m at a particular resolu- 
tion, then the largest £ retained in the right-hand side 
of the Ei equations after filtering is 2^ max /3 — 1, and 
the largest I retained in the right-hand sides of the 
equations is 2£ max /3. This filter is a variation on the 
standard 2/3 rule used to remove the inevitable effects of 
aliasing whenever functions are multiplied using spectral 
methods H3. 

For the first of these tests we choose initial data that 
corresponds to a static point charge that is located in the 
hole in the center of the computational domain. Thus 
we choose initial data with Ei — —di(f> — r~ 2 dir and 
— 0, appropriate for a unit point charge at rest at the 
origin. We then solve Eqs. I©-© with 71 = I/72 = 0.1 
on a computational domain with topology S 2 x R 1 de- 
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fined by 1.9 < r < 11.9. (This is the same compu- 
tational domain that we typically use to evolve single 
black hole spacetimes.) At both the inner and outer 
spherical boundaries we set the time derivatives of the 
incoming characteristic fields to zero, i.e. we impose 
freezing boundary conditions, Eq. I|38|) . The scalar po- 
tential 4> is held constant in time. We find that these 
numerical evolutions are stable and convergent and the 
constraints are preserved, as shown in Figs. and 0] 
These computations were performed with radial resolu- 
tions N r = 11, 21, 31, 41, 51, 61, and 73 collocation points, 
and a fixed angular resolution with spherical harmonic 
index £ max = 5 (or equivalently, Ng — 6 and N v = 12 
angular collocation points). For £ max = 9 the results are 
indistinguishable on the scale of Figs. |31 and 0] except at 
the highest radial resolutions, indicating that the radial 
and temporal truncation errors dominate, as expected 
for a solution with little angular structure. This is a case 
(as we shall see) in which a time-independent solution 
is not always the best test problem to investigate the 
constraint-preserving properties of a system of evolution 
equations. 

Finally, we examine a highly dynamical solution of the 
fat Maxwell system on the computational domain S 2 x R, 
defined by 1.9 < r < 11.9. For this solution we choose 
initial data with Ei = and 

A x = A y =A z =e-^-^ 2 / w \ (62) 

where rg = 6.5 and w — 1.0. The initial values of Dij are 
set to the numerically determined values of diAj. These 
initial data correspond to a pulse of radiation centered 
at r = tq. This pulse is neither spherically symmet- 
ric nor even axially symmetric, because the Cartesian 
components of the vector potential are set to spherically 
symmetric functions in Eq. (|62[1 : however, only a small 
number of spherical harmonics are sufficient to represent 
the full solution. The scalar potential is set to <j> = for 
these solutions, and we impose freezing boundary condi- 
tions, Eq. (|38|l . on the incoming characteristic fields. 

Figure [S] shows a convergence plot for this case, us- 
ing evolution parameters 71 = I/72 = 0.1 and i? max = 5, 
which confirms that the numerical solution is convergent. 
For £ max = 9 the results are indistinguishable on the scale 
of the figure. Figure [B] shows that significant constraint 
violations do exist in these solutions with seven differ- 
ent radial resolutions: the data from the three highest 
resolutions coincide on the scale of this diagram. Thus 
the constraints are violated, but the constraint energy is 
still convergent in these solutions. This indicates that 
the constraint violation is a property of the true solu- 
tion of the continuum equations with freezing boundary 
conditions, rather than an effect caused by a defective 
numerical method. The constraint violation appears to 
be generated as the initially constraint-satisfying waves 
interact with the boundaries, starting at about issi 

We have also computed evolutions for these dynamical 
initial data using negative values of the evolution param- 
eters: 71 = I/72 = —0.1; the results are depicted in 
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FIG. 5: Convergence test for a dynamical solution of fat 
Maxwell on S 2 x R using freezing boundary conditions and 
positive values of j a . Shown are norms of differences between 
solutions at different resolutions: solid curves use L 2 norms 
and dotted curves use L°° norms. 
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FIG. 6: Constraint violation for a dynamical solution of fat 
Maxwell on S 2 x R using freezing boundary conditions and 
positive values of y a . Shown is the constraint energy {£) 
divided by the norm of the derivatives of the fundamental 
variables. 

Figs.0and[H] Since the product 7172 is unchanged from 
the previous runs, the characteristic speeds of the system 
remain the same. And the definition of the constraint en- 
ergy £ (which depends on the ratio 71/72) also remains 
unchanged; so this allows us to meaningfully compare £ 
for the two cases. For j a < the fundamental evolution 
system, Eqs. © — ©> is strongly but no longer symmetric 
hyperbolic as in the 7 a > case. Figures and |H1 show 
that these evolutions with negative values of j a appear to 
be convergent, and have fractional constraint violations 
that are comparable with those in the positive 7 a case. 
However as illustrated in Fig. [51 the evolutions in the 
negative j a case have constraint violating instabilities. 
We note that we also find a numerical instability for 
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FIG. 7: Convergence test for fat Maxwell on S 2 x R us- 
ing freezing boundary conditions and negative values of j a . 
Shown are norms of differences between solutions at different 
resolutions: solid curves use L 2 norms and dotted curves use 
L°° norms. 



FIG. 9: Constraint violation for fat Maxwell on S 2 x R. Solid 
curves are for 7 a > while dotted curves are for "f a < 
evolution parameters. Seven different resolutions are depicted 
for each sign of y a , but only the lowest resolution curves are 
distinct at the scales shown. 
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FIG. 8: Constraint violation for fat Maxwell on S 2 x R us- 
ing freezing boundary conditions and negative values of "/ a . 
Shown is the ratio of the constraint energy to the norm of the 
derivatives of the dynamical fields. 

the 7a < case, and apparently for all cases in which the 
fundamental evolution system is strongly but not sym- 
metric hyperbolic. This numerical instability appears to 
be associated with the angular discretization. It grows 
exponentially in time, and becomes worse at higher angu- 
lar resolutions. However, for the angular resolutions we 
use, and for choices of "f a near the symmetric hyperbolic 
range, such as the case shown here, 71 = I/72 = —0.1, 
the numerical instability is negligible compared to the 
constraint- violating instability shown in Fig. Only by 
going to higher angular resolution can one see any non- 
convergent growth at all on the time scales we consider 
here. For £ max = 9 the instability is visible only at late 
times (t 200) for the highest radial resolutions in the 
quantities ||<$u||^ 2 and ||£it||!,oo, and is not visible in plots 



FIG. 10: Angular convergence test for a dynamical solution of 
fat Maxwell, showing numerical instability for 7 a < 0. Shown 
is a norm similar to ||<5it|| Ij 2 defined in Eq. I58H except that 
the differences are taken between quantities at two different 
angular resolutions and fixed N r = 73. Dotted curves show 
results for 71 = I/72 = —0.1, and are labeled by the two 
angular resolutions that are subtracted. Solid curves show 
that the same quantities for the case of 71 = I/72 = +0.1 are 
convergent. 



of {£). To construct a quantity that is sensitive to this 
instability, we repeated the runs shown in Figs. EHHl at 
angular resolutions ^ max = 5, 7, 9, and 11 and computed 
the norms of differences of the fundamental fields at ad- 
jacent angular resolutions. These norms are plotted in 
Fig. HOI 

We see no indication of this numerical instability for 
values of 7 a in which the fundamental evolution system 
is symmetric hyperbolic (for example, the solid curves in 
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Fig. ^]are convergent). For choices of j a very far from 
the symmetric hyperbolic range (such as 71 = I/72 = 
10), the instability grows much more rapidly and dom- 
inates the results. Although it is possible that the nu- 
merical instability can be cured by modifying our angu- 
lar filtering algorithm, for the purpose of this paper we 
simply consider only values of 7 a and angular resolutions 
for which the timescale of this instability is longer than 
that of other effects we wish to study. 



B. Active Constraint Control 

In this section we investigate the use of the active con- 
straint control methods described in Sec. IIIII For the 
case of the fat Maxwell system this active control con- 
sists of switching the sign of the evolution parameters 71 
and 72 to ensure that the constraint energy (£ ) does not 
increase. In the previous section we presented two sets of 
numerical evolutions without constraint control differing 
only by the signs of 71 and 72 . The characteristic speeds 
and the definition of the energy {£) were the same for 
these evolutions. Both evolutions were convergent on the 
time scale considered here (t < 100), and the fractional 
constraint violation was convergent and similar in these 
cases on the same time scale. We now investigate the 
possibility of switching between these two cases during a 
single evolution as a means of reducing the constraint vi- 
olations. The strategy is to monitor the quantity on the 
right side of Eq. (|35[) , and to change the signs of 71 and 
72 whenever necessary to keep the right side negative. 
This should ensure that the constraint energy norm (£) 
is always decreasing, so the constraints should remain 
satisfied. Note that this method should work only as 
long as our numerical solution satisfies the equation gov- 
erning the evolution of the constraint energy, Eq. (|35|l . 
Figure ^] illustrates for the case j a > that this equa- 
tion does remain satisfied to truncation error level for the 
runs discussed in Sec. IV Al the plot for j a < is similar. 
Consequently we expected good results from this active 
control method. 

Figure IT21 shows the constraint violation for a case in 
which the j a are allowed to change sign at every time step 
in order to control the constraints. The signs are changed 
only if the right side of Eq. (|35|l becomes positive, and the 
current value of (£ ) exceeds the value it had after the first 
timestep. The latter condition is intended to prevent sign 
changes that attempt to reduce the constraint violation 
to less than truncation error. Since this constraint con- 
trol method depends on Eq. (|35|l being satisfied, control 
should only be possible to truncation error level at best. 
Figure 1121 shows that the maximum value of the con- 
straint is smaller than for the uncontrolled case (plotted 
as a dotted line in Fig. ll2[) . However, the improvement is 
only an order of magnitude at best, even for the highest 
resolution run. Even more disturbing is the lack of con- 
vergence of the constraint norm. The fundamental fields 
do not converge very rapidly (if at all) either, as can be 
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FIG. 11: Violation of the constraint energy evolution equa- 
tion, Eq. 1351 , for fat Maxwell on S 2 x R with freezing bound- 
ary conditions and 7 a > 0. Plotted is the difference between 
{£) and the time integral of the right side of Eq. 1351 for each 
resolution. 
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FIG. 12: Constraint violation (solid curves) for fat Maxwell 
on S 2 x R with active constraint control at every time step. 
Dotted curve is the comparable uncontrolled case. 



seen from Fig.^J This lack of effective constraint control 
is confirmed in Fig. ^] which shows that Eq. l|35f) is not 
satisfied very well for this case. It appears that the active 
constraint control mechanism severely degrades the con- 
vergence of our numerical simulations in such a way that 
Eq. (|35|) no longer holds to the needed or expected accu- 
racy. Consequently the active constraint control method 
is able to reduce the constraint violations only by a small 
amount over the uncontrolled case. 

One effect that can destroy convergence in these tests 
is the fact that the control mechanism is applied indepen- 
dently for each resolution (as pointed out by Tiglio [23). 
Therefore, at a given value of t, the evolution equations 
used for one resolution can be different (because of the 
sign of 71 and 72) than the equations used for another res- 
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FIG. 13: Convergence test for fat Maxwell on S 2 x R, for 
active constraint control at every time step. Shown are norms 
of differences between solutions at different resolutions: solid 
curves use L 2 norms and dotted curves use L°° norms. 
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FIG. 14: Violation of the constraint energy evolution equa- 
tion, Eq. 1351 . for fat Maxwell on S 2 x R with with active 
constraint control at each time step. Plotted is the difference 
between {£) and the time integral of the right side of Eq. 135H 
for each resolution. 



olution at that time. This effect should have significant 
consequences only on quantities computed using differ- 
ent resolutions, such as the differences plotted in Fig. [Ill 
But quantities computed using a single resolution, such 
as (£), should not be affected. When these latter quan- 
tities are compared for different resolutions on the same 
plot, as in Figs. 1121 and 1141 the graph will not look like 
the "classic" convergence test in which all curves have 
the same shape. But the curves should (if everything 
else in the method is convergent) decrease at roughly the 
correct rate. Because we lose a great deal of accuracy in 
Figs . 1 1 21 and 1 1 41 compared to their uncontrolled counter- 
parts, Figs. Irjl and 1111 we believe that the fact that the 
control mechanism is applied independently for each res- 



olution is not the primary cause of the problem. In order 
to eliminate this effect on convergence, we repeated our 
simulations, but this time switched the sign of 7 a only 
once at the time t = 4 for each resolution, regardless of 
the sign of the right side of Eq. I|35[) or the magnitude of 
(£). In this case, exactly the same evolution equations 
are being solved at each resolution. As shown in Fig. 1151 
the convergence rate is severely reduced even in this case 
when the signs of -f a are switched at t = 4. Furthermore 
as shown in Fig. El Eq. i|35|) is violated after the signs 
are switched. 

We now believe that this nonconvergence is caused by a 
lack of smoothness of the fields that is introduced by the 
discontinuous change in the evolution equations: Sup- 
pose the signs of "f a are switched at a time to, and suppose 
that at time just before this switch (t = to — e) some out- 
going characteristic fields at the boundary are nonzero. 
When the signs of j a are switched, some of the outgoing 
and zero-speed characteristic fields will be converted to 
incoming characteristic fields, and vice versa, as can be 
seen from Eqs. (|13fl - l|17|) . (For example, switching the 
signs of the j a in Eq. El while keeping Ei and Dij fixed 
yields ?/f fter = — ^before-) Therefore, at a time just after 
this switch (t = to+e) the solution will contain some non- 
vanishing incoming characteristic fields near the bound- 
ary. However, our freezing boundary condition requires 
the incoming characteristic fields to be constant in time 
at the boundary. Since the incoming characteristic fields 
propagate inward, at times after the switch (t > to + e) 
a kink will appear in the profile of these incoming fields. 
This type of boundary-condition-induced kink is illus- 
trated in Fig. El The sketch on the left in Fig. E| illus- 
trates an incoming characteristic field at time just after 
the switch (t + e) , and the sketch on the right shows the 
kink in this field that develops from the boundary condi- 
tion. The existence of such a kink in the evolution fields 
greatly reduces the convergence rate of our spectral evo- 
lution method. And even for finite-difference methods 
such a kink is likely to reduce the convergence as well, 
since a kink in the fundamental quantities translates into 
a discontinuity in the constraints. Unless great care is 
taken to ensure that discontinuous solutions are treated 
properly (a standard problem in computational fluid dy- 
namics but quite foreign to vacuum numerical relativity 
because the gravitational field is not expected to have 
physical shocks), Eq. (|33|l is likely to be violated and the 
constraint preserving mechanism will fail. We have made 
several attempts to replace the freezing boundary condi- 
tion with a condition that smoothly adjusts the value of 
the incoming fields at the boundary. Unfortunately none 
of these attempts have been very successful. 



C. Constraint-Preserving Boundary Conditions 

Finally, we have performed a series of tests on the 
constraint preserving boundary conditions described in 
Sec. IIVI Figures El an d El show evolutions of our dy- 
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FIG. 15: Convergence test for fat Maxwell on S 2 x R, with 
the signs of *y a switched at t — 4 for all resolutions. Shown are 
norms of differences between solutions at different resolutions: 
solid curves use L 2 norms and dotted curves use L°° norms. 
Compare to runs with fixed y a in Fig. |7| 




Time 

FIG. 16: Violation of constraint energy evolution equation, 
Eq. 1351 , for cases with ~y a switched at t — 4 for all resolutions. 
Plotted is the difference between {£) and the time integral of 
the right side of Eq. 1351 for each resolution. Inset shows 
detail at early times, showing that Eq. 1351 is satisfied until 
the switch at t = 4. 

namical initial data on S 2 x R (analogous to that used 
in Figs. and IHJ) with the boundary condition on U 2 ~ 
now set according to the constraint preserving condition 
Eq. H50JI . The -f a are negative for the plots in Figs. 1181 
and 1191 The constraints are satisfied, and the simula- 
tion appears to be convergent (except for a late-time 
angular numerical instability, not visible on the plots, 
that appears identical to the numerical instability dis- 
cussed at the end of Section IV All . We have also per- 
formed these evolutions using positive values of j a , and 
the simulations appear to be convergent, completely sta- 
ble, and constraint-preserving in this case. Figure [201 
compares the unnormalized constraint energy for these 




FIG. 17: Left curve represents a characteristic field at one 
instant of time, and right curve the evolution of this field 
at a later time. Freezing boundary conditions produce the 
non-smooth but continuous solid curve extension, while stan- 
dard maximally dissipative boundary conditions produce the 
discontinuous dashed curve extension. 



evolutions with those run with freezing boundary con- 
ditions. Thus adopting constraint preserving boundary 
conditions clearly does improve the constraint preserving 
properties of these evolutions much more than the active 
constraint control method. 

But this is not the entire story. Figure [2] shows the 
norm of the fundamental dynamical fields ||w|| 2 for evo- 
lutions using constraint preserving boundary conditions 
with 71 = —0.1 (solid curves) and 71 = 0.1 (dotted 
curves). This plot shows that while the positive 71 evolu- 
tions are stable, those with negative 71 are not. A more 
extensive sampling of the parameter space reveals that 
evolutions preformed with 7! = I/72 = {0.1, 1.0, 2.5, 2.9} 
(for which the principal evolution system is symmetric 
hyperbolic) appear to all be convergent, constraint pre- 
serving, and stable. Conversely, we find that evolutions 
performed with 71 = I/72 = { — 1.0,-0.1,3.5,4.1} (for 
which the principal evolution system is strongly but not 
symmetric hyperbolic) are all constraint preserving but 
unstable. These evolutions are numerically convergent 
for the resolutions and time scales we have tested (ex- 
cept for a slowly-growing angular numerical instability, 
not visible on Figure |^ that appears at late times or for 
high angular resolutions). Therefore, the type of growth 
seen in Figure |2] appears to represent a solution of the 
partial differential equations. Since these solutions do 
satisfy the constraints, the driving force for these in- 
stabilities must be an excess of incoming radiation that 
is reflected back into the computational domain by the 
boundary condition. We refer to this type of instabil- 
ity as a boundary condition driven instability. Thus the 
constraint preserving boundary conditions are a great im- 
provement over the other methods studied here, but they 
do not completely eliminate all the instabilities in these 
strongly hyperbolic cases. Further study will be needed 
to determine whether these boundary conditions can be 
improved. 
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FIG. 18: Convergence test for fat Maxwell on S 2 x R, us- 
ing constraint-preserving boundary conditions and y a < 0. 
Shown are norms of differences between solutions at different 
resolutions: solid curves use L 2 norms and dotted curves use 
L°° norms. 



FIG. 20: Constraint violation for fat Maxwell on S 
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FIG. 19: Constraint violation for fat Maxwell on S 2 x R us- 
ing constraint-preserving boundary conditions and 7 a < 0. 
Shown is the constraint energy {£) divided by the norm of 
the derivatives of the fundamental variables. 



VI. DISCUSSION 

This paper explores the effectiveness of two methods 
for controlling the growth of constraints in hyperbolic 
evolution systems. Using an expanded version of the 
Maxwell evolution system — which we call the fat Maxwell 
system — we showed that significant constraint violations 
and in some cases even constraint violating instabilities 
occur when the evolutions are performed using "stan- 
dard" numerical methods and boundary conditions. We 
show that the active constraint control method (which 



d(wl 

m 



has been studied by Tiglio and his collaborators [23,123]) 
is not very effective in controlling the growth of con- 
straints in the fat Maxwell system when spectral numer- 



for 7 a < 0. Solid curves use constraint-preserving boundary 
conditions while dotted curves (same as the dotted curves in 
Ffg-EJ use freezing boundary conditions. Seven different res- 
olutions are depicted for each type of boundary condition, 
but only the lowest resolution curves are distinct at the scales 
shown. 
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FIG. 21: Norm of the fundamental variables for fat Maxwell 
on S 2 x R with j a < (solid) and "f a > (dotted), using 
constraint-preserving boundary conditions. Even though the 
constraints are satisfied for f a < 0, the fundamental quanti- 
ties increase exponentially, but in a convergent manner. Seven 
different resolutions are depicted for each case, but only the 
lowest resolution curves are distinct at the scales shown. 



ical methods are used. This lack of effectiveness appears 
to be caused by the non-smooth nature of the control 
mechanism for this system. We also show that constraint 
preserving boundary conditions are very effective in sup- 
pressing the constraint violations in this system. Un- 
fortunately these constraint preserving boundary condi- 
tions did not eliminate the instabilities for the strongly 
(but not symmetric) hyperbolic evolution systems. In 
these cases these boundary conditions merely converted 
a constraint violating instability into a boundary condi- 



15 



tion driven instability. Generalizing these methods to 
more complicated systems like the Einstein evolution 
equations should be straightforward. For more general 
systems the analogue of the constraint energy evolution 
equation, Eq. 150|) . will contain both boundary terms like 
the fat Maxwell system and also volume terms, so that 
constraint violations can be generated both at bound- 
aries and in the bulk. We expect that constraint pre- 
serving boundary conditions alone will not be sufficient 
to control the constraint violating instabilities that occur 
in these more general systems. Instead we expect that 
some combination of methods will be needed. The disap- 
pointing results obtained here with the active constraint 
control mechanism suggest that significant improvements 



will be needed to make this method useful for helping to 
control the growth of constraints in the Einstein system 
for evolutions based on spectral methods. 
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